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Abstract. By using the numerically exact density-matrix renormalization group (DMRG) 
approach, we investigate the ground states of harmonically trapped one-dimensional (ID) 
fermions with population imbalance and find that the Larkin-Ovchinnikov (LO) state, which 
is a condensed state of fermion pairs with nonzero center-of-mass momentum, is realized for 
a wide range of parameters. The phase diagram comprising the two phases of i) an LO state at 
the trap center and a balanced condensate at the periphery and ii) an LO state at the trap center 
and a pure majority component at the periphery, is obtained. The reduced two-body density 
matrix indicates that most of the minority atoms contribute to the LO-type quasi-condensate. 
With the time-dependent DMRG, we also investigate the real-time dynamics of a system of 
ID fermions in response to a spin-flip excitation. 
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1. Introduction 



One-dimensional atomic gases have been vigorously investigated over the last decade. 
By using a strongly anisotopic Ioffe-Pritchard-type magnetic trap HI or pairs of 
counterpropagating laser beams [2], elongated traps of atoms have been realized. Tens 
to hundreds of ultracold atoms have been trapped in such an elongated potential that the 
maximum kinetic energy of the atoms are much smaller than the energy separation between 
the two lowest transverse levels of the potential J3l SI |5l [6l 13. Thus the physics of the 
system becomes essentially one-dimensional [[81 |9]|. For the case of bosons, a strongly 
interacting Tonks-Girardeau regime was realized [3J; no equilibration of an integrable system 
— "quantum Newton's cradle" — was demonstrated 0; and the crossover from thermal to 
quantum noise was also observed [6]. 

For the case of fermions, the ID atomic gas allows one to study an analogue of a system 
of electrons in solids. The numbers of atoms in trappable hyperfine states can be controlled 
individually, so that by populating two hyperfine states simultaneously, one can study effects 
of population imbalance on the properties of one-dimensional systems. While a true long- 
range order is forbidden in ID even at T = by the Hohenberg-Mermin-Wagner theorem, 
(quasi) condensation and superfluidity are possible in a finite-sized, trapped system. The 
search for Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) IfTOlfTTTl pairing states has been of major 
experimental interest in ID fermionic systems. 

The confinement-induced two-particle bound state was observed flU in an array of quasi 
one-dimensional tubes with equal populations of two hyperfine states, validating the ID 
treatment of these tubes. More recently, the Rice group has realized quasi one-dimensional 
systems of population-imbalanced cold atoms and measured the density profile for each of the 
hyperfine species. 

In this article we first report our study on the ground state of a system of population- 
imbalanced fermions in harmonic traps [12]. The system has also been numerically studied 
by using the density-matrix renormalization group (DMRG) lfT3l fl4l [T51 IT6l [T71 IT8l [T9ll. 
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Bogoliubov-de Gennes (BdG) equations ESllIIl, and Quantum Monte Carlo (QMC) method 
E2l 1231 |24| . The phase diagram including an FFLO superconductor region for a ID 
homogeneous system was obtained by bosonization ll25l . Studies on trapped systems based 
on the exact solution of the Gaudin-Yang model were carried out in Refs. [|26ll27l[28ll29l[30ll . 
Other related studies include the cases with unequal (effective) masses and the balanced or 
imbalanced numbers of atoms ll3Tl[32l[33l[34ll : an angular FFLO state in a toroidal trap 11351 : 
effects of particle-correlated hopping in an optical lattice ll36lk imbalanced Fermi gases on 
two-leg ladders OTll and two-dimensional optical lattices 113811 : composite particles of more 
than two particles [f39ll : universality and Efimov states in systems where two or more species 
of atoms are confined in traps having different dimensions ll40l |4T| : transport through a Y- 
junction and ring-geometry systems Il42ll . 

We note that, for the ground state that is symmetric under time reversal, the condensate 
cannot be of a Fulde-Ferrell type, because the Fulde-Ferrell state is described by a pair 
correlation function that is a complex- valued function of the spatial distance. On the other 
hand, the Larkin-Ovchinnikov state is not forbidden by the time-reversal symmetry because 
the correlation functions can be real-valued. 

Here we numerically study the ground-state wavefunction of the trapped system by 
discretising the system and the DMRG R3l 1441 1451 method. We calculate the density 
distribution of individual spin components in the trap, which is a physical quantity measured 
in in-situ light absorption or phase shift measurements Il46ll47ll48l . and show that the density- 
difference distribution exhibits oscillations which are the hallmark of the Larkin-Ovchinnikov 
state Il49ll . 

In the second part of the paper, we study the real-time evolution of a system of trapped 
Fermi atoms after a sudden spin flip. Because of the large interparticle distance and low speed 
of the atoms, the dynamics is extremely slow as compared to electrons in solid state physics. 
Moreover, the atoms are trapped in an optical and/or magneto-optical potential in the vacuum, 
so that the dynamics can be optically observed in real time. As the state of the atoms can be 
manipulated by microwaves, it is of interest to study the response of the ID ground state to a 
change in population. 



2. The ground state 

First, we study the ground state of the trapped, population imbalanced fermionic system in ID 
within the single-channel model. The effect of a molecular state in the BCS-BEC crossover 
in ID systems ll50l has also been studied in population imbalanced systems [[18115111 . but here 
we focus on the case in which the contribution of the molecular state can be neglected. 

The interaction between ultracold atoms is described by the s-wave scattering which can 
be approximated by the Lee-Huang- Yang pseudopotential [l52l 

U{r) = gf5{f)hr-), (1) 
or 

with the bare 3D coupling constant g^ D . The relation between gl° and the s-wave scattering 
length af* is given by gf = InTiarf /pi, where /a is the reduced mass. 
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The relation between a 3 ® and the scattering length in a ID trap a\ D has been given in S, 
for the case in which the transverse confinement can be approximated by a 2D harmonic 
potential with (angular) frequency a> ± . Provided that the size of the ground state of the 
transverse confinement a ± = (n/fuo ± ) 1/2 is much larger than an effective range r of the bare 
potential, a' D is given by 



a 2 I a 3D 

ID _ u ± ' - ~ 



with C = 1.4603 The ID scattering amplitude for the relative wavenumber k z is expressed 

as 

feven(k z ) = — - ■ - TZTTj ToT* (3) 

1 + ik z a; u + O {(k z a ± Y) 

which, in the low-energy limit (\k z a ± \ «: 1), reduces to the scattering amplitude derived from 
the one-dimensional 8 potential 

U m (z) = gl B S(z) (4) 

withgJ D = -h 2 l{puxf)M. 

When we have two fermionic species with the equal mass m , the reduced mass is 
fji = m /2. The two hyperfine states can be regarded as a pseudo spin- 1/2 system with spin up 
(| T)) and down (| D) states. We assume that the up-spin state is the majority state. We start 
with an effective ID Hamiltonian for the continuum, 



n-f = 



n rr > a ->V ' n n' 



where m is the atomic mass, Viz) = kz 2 /2 is the harmonic potential in the z direction and g is 
the coupling constant. We discretize this Hamiltonian to perform calculations in the DMRG 
framework. 

We take a characteristic size of the potential, /, as the unit of length, and denote a 
characteristic trap depth by A: V(+l) = A. We discretize the system by introducing a lattice of 
L sites with lattice constant d = 21/ L. The transfer amplitude between the neighboring sites 
is chosen to be J = ti 2 /(2md 2 ). The band dispersion, E(k) = 27(1 - coskd) = J(kd) 2 + 0(k 4 ), 
where k = p/fi is the wave number of the atom with momentum p, reproduces the energy 
dispersion of the atom in the free space, E(k) = p 2 /(2m) = Ji 2 k 2 /(2m), in the vanishing limit 
of the filling factor (L — » oo). The interaction g iD S(za,-[ - Z a ',\) is approximated by an on-site 
interaction, U 2^jo "i,t"U> wnere z runs over ai l lattice sites 0,1,..., (L— 1), with coupling 
constant U = g iD /d. In the following we take h = m = 1. We consider the case of negative 
U, which corresponds to a positive a| D and attractive interaction between atoms with opposite 
spins. 



2.1. Relations between the experimental and discretized model parameters 

Here, we describe how the various parameters introduced above can be related to the 
experimental parameters. For the 40 K atoms with a> r = 2n x 69 kHz and a> z = 2n x 256 Hz 
flU, a ± = 86 nm, and at a 3D Feshbach resonance (a^ D — » oo), we have a] D = 62.8 nm. 



A typical size of the axial ground-state wave function of a single atom in this trap is a z = 
(h/nu) z y /2 = 1.41 yum. As the (half-)width of the simulated region of the harmonic trap, we 
take / = 6.28 /mi. The length will be measured in units of /below. Now, 7 = h 2 / (2m(2//L) 2 ) = 
h 2 L 2 /(Sm f), so that A/7 = (8Am / 2 lh 2 )L~ 2 and U/J = -n 2 /(m a] D d/2)/J = -(8//a 1D )/L. 
We choose A = k B x 246 nK, so that A/7 = 6400L" 2 and U/J = -800L" 1 , unless otherwise 
stated. 

If there are N atoms per spin, the Fermi energy for the noninteracting case is E F = Nha> z . 
For the same number of atoms per spin, the dimensionless strength of interaction £ at the 
Fermi level is 



m = (k F a?y l = vw^^ = U2N - 



U 1 - 



(6) 



Similarly, for n atoms per spin per unit length /, the Fermi wave vector is k P = nn/l, and we 
obtain 

£(n) = (k F a]°y l = (l/nal D )n~ l = 31.8m" 1 . (7) 



2.2. Method: Density-matrix renormalization group 

After the discretization described above, the Hamiltonian is given by 

L-l L-l L-l 

cr (=1 i=0 !=0 

where V(i) = k(J - C) 2 with k = 4A/L 2 and C = (L - l)/2. Note that A/7 cc L' 2 (thus 
/c/7 oc LT A ) and U/Joc L~ l . We use DMRG to calculate the ground state of the Hamiltonian 
([8]) within the particle-number sector of and Ni atoms in spin-up and spin-down states, 
respectively. The population imbalance parameter is defined as P = (iVj - Ni)/N tota i, where 
Motai = + TYj,. The Fermi momentum at the trap center is calculated from the averaged 
density n^ = {h i tT ld) as k Fcr = rm^, where the overbar denotes the average over 0.1L - 0.2L 
neighboring sites. 

Typically m = 300 states per block are retained and 4-5 finite system iterations are 
required for DMRG convergence after the standard infinite system initialization or the 
recursive sweep method initialization [|53l . We target the ground state with 7V T and N± atoms 
in spin-up and spin-down states as well as several states with a few atoms added or removed, 
in order to improve the convergence of the pair correlation. While we do not assume the 
reflection symmetry of the wavefunction in the finite system iteration process, the reflection 
symmetry emerges spontaneously for the ground state wavefunction for any of the set of 
parameters studied. The sum of the eigenvalues of the reduced density matrix for the discarded 
eigenstates in each DMRG step is below 2 x 10~ 5 for N = 64 and L = 320; it is below 5 x 10~ 7 
for N = 40 and L = 200 if we target the ground state alone. 



6 



2.3. Results: the Larkin-Ovchinnikov quasi-condensate 

2.3.1. Pair correlation function We first examine the on-site pair correlation function 
defined as 



where O on - s i te (zj,Zj) = {n^n^) gives the double occupancy ratio. The left column of figured] 
shows O on _ s i te (z 7 -, zj) in color-coded two-dimensional plots for several different values of the 
imbalance parameter P. For P = 0, the pair correlation is maximal for z, = Zj with the double 
occupancy per site, decreases with increasing \zi - Zj\, and precipitously drops to zero where 
the atomic density, shown in the right column, vanishes. 

As we increase P, the correlation function starts to oscillate. The zeros of the pair 
correlation O n~site(z> zj) for fixed zj appear almost periodically, except for the region \z-Zj\ <^ 
/, where O on -site is positive and large. The sign change of the pair correlation reflects that of 
the pair amplitude. We note that the whole region with a significant density of minority atoms 
is strongly correlated inside the ID trap. 

2.3.2. Density distribution: Oscillation of density difference 

Next, we examine the density distributions of the spin-up and spin-down atoms. For 
the non-interacting system, the extensions of the majority and minority distributions are 
determined by those of the A^-th and A^-th eigenstates of the trap, respectively. As \U\ is 
increased, the attractive interaction causes the atomic distributions to shrink towards the trap 
center, as shown in figure |2] We also find that the density difference (n^ - n^/d is initially 
spread all over the region with a nonzero atomic distribution, but that it shrinks more rapidly 
than the atomic distribution for sufficiently large \U\. The oscillation is independent of the 
lattice constant d for large L, as can be seen in figure |3] The peaks in the frequency spectrum 
D(q) of the density difference, calculated as 



and shown in the right column of figure [31 do not show a significant shift upon increasing the 
system size. 

From the spin-dependent Fermi momentum at the trap center, k Fcr , we obtain the 
difference of the Fermi momenta, q = k^ - kpi- In the LO state, an up-spin atom close to the 
up-spin Fermi point (in ID) +kp\ and a down-spin atom close to the down-spin Fermi point 
+&fi form a pair; therefore the total momentum of the pair must be close to +q. In figure |4] 
we plot the cosine transformations of the pair correlation O on - s ite(z, Zl/z), where zl/z = l/L 
is the location of one of the two sites closest to the trap center z = 0, and the density 
difference 2S z (z) = n^(z) - n±(z) against the wave vector. To clearly identify oscillations close 
to the trap center, the density difference has been multiplied by a slowly varying numerical 
factor, g(z) = exp (-6z 2 ), before Fourier transformation. The Fourier components of the pair 
correlation and density difference at wave vector q are defined as 



O, 



on 



site(Zi,Zj) = (^olc^C^Cj^C jil^o), 



(9) 




(10) 




(11) 
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Figure 1. Left column: pair correlation function © for L — 200, N = 40, A/ J = 0.16 = 
6400/L 2 , U/J = -4 = -800/L, and P = 0,0.1,0.3,0.7. Right column: the corresponding 
density distribution functions for majority (up), minority (down), and their difference (diff.) 
with (int.) and without (nonint.) interaction. (a)(b): reproduced from |[T2l with modified color 
codes. 



and 



n<tis(q) 



2cos(fei)f(z)2S z fe) 



(12) 



As shown in figure HI the oscillation wave vector of the pair correlation function is 
always close to q. This is consistent with studies of population imbalance in the optical 
lattice systems, namely, the pair momentum distribution function calculated in |[T3Tl . Fourier 
transform of the pairing correlation Ifl4l . on-site pair amplitude lfT5l . spatial noise correlations 
[fT6l , and radio-frequency spectroscopy calculation Il2~0~ll , as well as the pair momentum 
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Figure 3. Left column: Density distributions of majority and minority atoms and their 
difference plotted against the normalized coordinate in the harmonically trapped system of 
fermions with (Nf,N±) = (24, 16) atoms and for different values of L. A/ J = 6400/L 2 and 
U/J = -800/L. Right column: Fourier (cosine) transformations of the density distribution 
plotted against the wave number. 



distribution calculated in Il22l |23l for the continuum. The oscillation wave vector of the 
density difference is close to 2q, as expected in the LO state ||49T| . This result is consistent 
with the spin-spin correlation function studied in Il24l and distribution of spin polarization 
studied by solving BdG equations in [|2T|. 



2.3.3. Phase diagram: Parameter region with phase separation We find the following two 
phases: i) the spin-down (minority) atoms occupy a smaller region compared with the spin- 
up (majority) atoms and the periphery of the system is occupied by the spin-up normal 
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Figure 4. Upper : Fourier components of the pair correlation function (left) and the density 
difference (right; see the main text) plotted against the wave vector. Lower (reproduced from 
ifPUl with modified color codes): The wave numbers of the oscillations of the pair correlation 
function and the density difference plotted against the difference of the Fermi wave numbers 
at the trap center. The values of the parameters are L = 200, N = 40, A/ J = 0.16 = 6400/L 2 , 
U/J = -4 = 800/L. 



component alone, and ii) the densities of spin-up and spin-down atoms are different at the 
trap center but become almost equal over a finite width close to the periphery. These two 
situations can be thought of as two "phases" of the system. The transition between the two 
phases occurs when the trap is filled with the LO condensate. Such a phase transition was 
predicted in studies in which the exact solution of the Gaudin-Yang model without a trapping 
potential was plugged into the trapped system via the local-density approximation (LDA) 

For a given scattering length a or coupling constant g, case i) is realized for greater 
spin polarization P, where the condensate can no longer hold all the majority atoms and the 
system gains energy by pushing some of the majority component to the edge of the trap. This 
is because the condensation energy of the remaining LO condensate increases by a decrease 
in the Fermi wave-number difference and the number of nodes of the pairing potential. 

Because the coupling constant \g\ increases with increasing \U\, which corresponds to 
shorter \a\ in ID, the critical point of phase separation shifts toward greater values of P, as 
shown in figure |5J The phase separation curve, however, does not exceed P 0.21. 
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Figure 5. Phase diagram of the trapped population-imbalanced Fermi gases: i) an LO 
condensate at the trap center and spin-polarized normal gases that include only the majority 
component at the edges of the trap, and ii) an LO condensate at the trap center and equal- 
population condensates at the edges of the trap. 
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Figure 6. Density distributions for spin-up (top) and spin-down (bottom) atoms plotted against 
the normalized trap coordinate shifted by zo, where zo is the solution of V(zo) = 0, so that the 
z - Zq — corresponds to the trap minimum. The parameters used are L — 200, N = 40, 
P = 0.2, A/ J = 6400/L 2 = 0.16, and U/J = -500/L = -2.5. 



2.3.4. Effect of asymmetric potential: Accuracy of local density approximation 

To study the robustness of the Larkin-Ovchinnikov state with respect to the trap geometry 
and the slope of the trap, we examine the ground state of the system for a spatially asymmetric 
potential, 

asymU,j ~ \ kM-C) 2 (i<C) ' K S) 

where 0<C<L-l,/c + = A/C 2 and k~ = A/(L - 1 - C) 2 . The bottom of the trap is located 
at i = C = L V?/( ^ + 1). Let r = k + /k- be the ratio of the spring constants of the harmonic 
trap; this ratio equals the ratio between the potential slopes in the right and left sides of the 
bottom of the trap. Figure [6] shows the density distribution against the shifted trap coordinate 
z - Zo, where zo = 2C/L- 1 = (V^ _ 1)/(V^ + 1) and V(zo) = 0. Figure [7] shows the same 
density distribution plotted against the site energy level V{z). 

Despite the increased asymmetry of the trap, the atomic distribution changes very little, 
and the pair correlation exhibits similar oscillations as those shown in figure [8] This is yet 
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Figure 7. Density distributions for spin-up (top) and spin-down (bottom) atoms for L = 200, 
N = 40, P = 0.2, A/J = 6400/L 2 = 0.16, and U/J - 500/L = -2.5, plotted against the site 
energy level V(z). Here r = k+/k- gives the ratio between the potential slopes shown in dT3b . 




Figure 8. On-site pair correlation for r = 1,4,9, 16, L = 200, N = 40, P = 0.2, 
A/J = 6400/L 2 = 0.16, and U/J = -500/L = -2.5 shown in a color-coded plot against 
the normalized trap coordinates. 

another evidence that the LDA is a rather good approximation in a system of ID population- 
imbalanced fermions. Recently the phase separation between the FFLO condensate and fully 
paired or fully polarized wings of the trap has been studied in |fT9l , where a completely 
asymmetric setup, which corresponds to r — » co, was used and compared with analytic results 
based on the combination of the Bethe ansatz solution and LDA. 
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Figure 9. Largest eigenvalue A™ x of the reduced two-body density matrix plotted against the 
imbalance parameter P for N = 20, 30 and 40. The parameters are L = 80, A/ J = 1 = 6400/L 2 
and U/J = -10 = -800/L. Reproduced from lfT2l with modified color codes. 

2.3.5. Two-body pair correlation matrix 

To quantify how much atoms contribute to the (quasi-) condensation in the population- 
imbalanced ID system, we examine the behaviour of the eigenvalues of the two-body reduced 
density matrix, from which we can identify the presence or absence of the off-diagonal long- 
range order [l54ll . We consider the two-body reduced density matrix 

pffx* = (^plc^o), (14) 

where a,fi, y, 8 are one of the M possible states that N fermions in the system can each take. 
As proved in 11541 . the upper bound of the eigenvalue A 2 of p f ^. yS is N(l —(N — 2)/M), and the 
system possesses an off-diagonal long-range order if and only if A 2 is of the order of N. 

Since pairing occurs between up-spin and down-spin atoms, it suffices to calculate the 
part of p Ml for which a and f3 as well as y and 6 have different spins. Therefore, we calculate 

Pi h h;j r ,jl = (^olc^C.^Cj^Cj^o). (15) 

The eigenstates of this matrix with large eigenvalues correspond to states with high pair 
concentration. It can be shown that they become independent of the lattice discretization 
in the limit of L — » oo. Note that while the dimension of the matrix is L 2 , the trace of the 
matrix, which is also the sum of the eigenvalues, equals A^A^ <sc L 2 , and the maximum 
possible eigenvalue, which is close to N±, is realized when every spin-down atom is paired 
with a spin-up atom and when these pairs are all condensed into a single state. 

We have diagonalized the two-body density matrix for the ground state that is obtained 
in the DMRG simulation for each parameter set. Because the number of independent matrix 
elements scale as 0(L 4 ) with increasing the number of sites L, we used L = 80 and limited 
N up to 40. The number of states per DMRG block is also limited to m = 100, because for 
L = 80, the energy and largest eigenvalues of the two-body density matrix for the ground state 
converges for smaller m compared to the L = 200 case. The i = i', j = f components of this 
matrix are the on-site pair correlation functions and correspond to the pair amplitude of the 
condensation of on-site "molecules." In Ref. [13], the matrix of the on-site pair correlation, 

P™ = Oon-sitefo, Zj) = <%|cJ/J/;, T C iU |¥ >, (16) 

was diagonalized and the resulting eigenfunctions showed periodic oscillations. A similar 
feature is confirmed in our analysis of the full two-body pair correlation matrix. 
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Figure 10. Left column: Distribution of the «-th largest eigenvalues of the two-body density 
matrix p <fl~5T > plotted against P for N = 40. The observed kinks in the distribution curves 
are marked with vertical arrows. The parameters are L = 80, A/J — 1 = 6400/L 2 and 
U/J = -10 = -800/L. Right column: Numbers of majority and minority atoms, and the 
sum of the eigenvalues of p that are enhanced by the effect of pair quasi-condenstaion, plotted 
against P, for L = 80, N = 40, A/J = 1 = 6400/L 2 and U/J = -10 = -800/L. 



In figure we plot the largest eigenvalue of p against the imbalance parameter P for 
different numbers of atoms N. The dependences on P for different values of /V are strikingly 
similar. This indicates that, in the trapped system, the lowest energy state of atom pairs is 
occupied only by a few pairs, before /V reaches 20, and they do not accommodate additional 
pairs of atoms, probably due to the effect of the trap potential which localizes the low energy 
states. 

In the left column of figure [TO] we plot the n-th largest eigenvalue of p against n. While 
the distribution is by definition a decreasing function of n, we find a kink in the distribution 
at different values of n for different P. Without the interaction, the eigenvalues of the reduced 
density matrix are always unity or zero, because each single-atom state is either filled or 
empty. The kink in the distribution suggests that the eigenvalues larger than the kink value 
are enhanced by the effect of pair quasi-condensation. In the right column of figure \\0\ we 
plot the sum of such enhanced eigenvalues against P. The sum is close to the number of 
minority atoms in the system, which indicates that almost all the minority atoms contribute to 
the quasi-condensation and the oscillating pair correlation. 

3. Dynamics 

Now we study the dynamics of a ID, harmonically trapped Fermi gas after a population 
imbalance is introduced by flipping the spin of a single atom. The flipping can be realized by 
applying a microwave radiation, whose energy matches the hyperfine energy level separation. 
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3.1. Formulation: Real-time DMRG simulation after a spin-flip excitation 

We simulate a sudden spin flip by applying a linear combination of local spin raising operators 



under spatial discretization. We obtain the wavefunction of the system for the ground state in 
the Nf = Ni = Af-spin sector |T ) as discussed in the previous section, and then apply S + to 
obtain the normalized (with a constant AO wavefunction as the initial state of the spin-flipped 
system, 



The spin flip does not in general map the ground state in the (A^, N\) = (N, N) sector to 
an eigenstate in the (A^ ,N±) = (N +1,N -1) sector. Therefore, the system will evolve in time, 
showing some real-time dynamics. As a realistic case, we consider a situation in which the 
coefficient c/ is spatially localized, which corresponds to the local spin excitation. This can be 
implemented experimentally by using a magnetic field gradient. 

In the case of an attractively interacting system, the spin flip leads to pair-breaking and 
increases the energy of the system relative to the ground state in the (N+ 1, N- 1). As discussed 
above, this ground state consists of a Larkin-Ovchinnikov condensate at the trap center and 
an equal-population condensate near the trap fringes, due to the low (P = l/N < P c ) spin 
polarization. 

For the exact calculation (after the spatial discretization), the full Hilbert space has to 
be considered, and the size of the Hilbert space rapidly increases. For example, for N = 4 
in the L = 16-site system, the dimension of the Hilbert space is already 2446080. Several 
methods to apply DMRG to simulate real-time dynamics of quantum systems have been 
implemented j53 [56l 1571 l58l |59l l60l l6Tl [62l . Comparison of some of these methods have 
been conducted in [|63l . From the ground state in the restricted Hilbert space obtained in the 
DMRG simulation, we calculate the spin-flipped state, and take this state as the state at t = 0. 
Then we apply the time-stepping targeting method lloTI to simulate the real-time evolution of 
this initial wavefunction \¥(t = 0)). 

In the time-stepping targeting method ll6TTl . the Schrodinger equation 



is treated as the equation of motion for the many-body wavefunction |*F(0), and the Runge- 
Kutta method is used to obtain |*F(? + At)), where At is the time step, within the Hilbert space 
selected in the previous finite-size system DMRG step. At each step, |*F(f)), + At)), 
and several other wavefunctions ^(t + ajAt)) are targeted in the choice of basis. After 
one iteration of the finite-system algorithm DMRG, t is increased by At and the process is 
repeated. Following flU, here we target four time steps |¥(0>, PP(* + Af/3)>, \*¥(t + 2Af/3)>, 
and |*P(? + At)), whose weights in the density matrix are 1/3, 1/6, 1/6 and 1/3, respectively. 
In this method, as in the methods based on the Suzuki-Trotter decomposition 11571 [58l . the 
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Figure 11. Energy of the simulated wave function after the spin-flip excitation at t — for 
L - 32, z c = 0, w = 1/16, U/J = -8 and A/ J = 1. E (N h Ni) denotes the ground-state energy 
in the sector of iVf spin-up atoms and Ni spin-down atoms of the Hamiltonian dS). The curve 
for Af = 0.001 is replotted in the inset. 



reduced Hilbert space in the subsystems after several DMRG iterations no longer contains 
much of the information on |*F(? = 0)). Use of the Runge-Kutta method in simulating real- 
time dynamics has also been discussed in II571 . 

Here we note that the dynamics of the superfluid to Mott insulator transition ll64l 1651 . 
spin-charge separation of a ID Fermi gas ll66l l67l . dynamics of spinless fermions after an 
interaction quench Il68ll69~l . expansion of a (bosonic or fermionic) Mott insulator [170117111721 . 
relaxation of atoms after a superlattice is modified Il73l 1731 . quench dynamics of the Bose- 
Hubbard model 11751 . and damping of dipole oscillations of a ID Bose gas after a sudden 
displacement of the trap 117611 have been studied by the time-dependent DMRG or a related 
method, time-evolving block decimation (TEBD) [1771 . For a mixed state, the real-time 
evolution of the density matrix of the system can be simulated in DMRG. Simulations of 
the dynamics for finite-temperature spin systems have also been carried out [17811791 . Here, 
we show the results of our simulation for the dynamics of pure states at T = 0. 



3.2. Results: Diffusion of the excess majority atom and destroyed correlation 

In the following, a spin-flipping excitation is parameterized by two parameters, the center of 
excitation z c and width w. The excitation factor for each site /, s h is determined by 

s, = Ccxp[-(z,-z c ) 2 /w 2 ], (20) 

where C is the normalization constant. 

First, we examine the case in which the spin is flipped close to the trap center (z c = 0). 
We consider the initial state obtained by flipping one of the down spins from the ground state 
in the (iV t ,iV 4 ) = (8, 8) sector, for U = -8 and A = 1. We take L = 32 and w = 1/16 and 
retain m = 300 states per DMRG block. Here, J = n 2 L 2 /(Sm l 2 ) = 39.4 nK x k B , and we take 
h/J= 194 //s as the unit of time. In figure[TT]we plot the evolution of the energy, defined by 

e(o = mowmt)), (2i) 
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Figure 12. Density distributions for different time steps after the spin-flip excitation at t — 
and the trap potential plotted against the normalized coordinate, for L = 32, z c = 0, w = 1/16, 
C//J = -8 and A/ J = V(+l)/J = 1. 

where is the calculated wave function at time f after the spin flip at t = 0. We observe 
that the energy goes up rapidly for At = 0.01 after t 1.3. Because the Hamiltonian is always 
the same, this increase in energy is due to the accumulated error in the Runge-Kutta numerical 
integration. For At = 0.002, the increase of the energy occurs later and it is less significant. 
For At = 0.001, the energy is almost constant until t 4. We present the result for At = 0.001 
up to t = 4. 

For U = -8, we plot the density distribution for different time steps in figure[[2l Because 
of the spin flip, initially (at t = 0) the majority density is largest close to the trap center, and the 
minority density has a deep dip there, making the density distribution localized there. As the 
time elapses, the majority density decreases and the minority density increases at the center, 
and the density difference exhibits two peaks at both sides of the central, decreasing peak. Up 
to t = 4, these two peaks move away from the trap center. There remain three peaks of the 
density difference. 

This is to be compared with the density distribution for the (Nf, N^) = (9, 7) sector of the 
Hamiltonian © (see figure [TBI. The majority and minority density profiles are not as smooth 
as those after the spin flip. Such oscillations in the density profiles have also been observed 
for larger L; they originate from the formation of tightly bound pairs, rather than the lattice 
discretization 031. In the excited state after the spin flip, such oscillations in the ground state 
are smeared out because of the weaker binding between spins and the excess kinetic energy. 

In figure \\3\ we see two peaks of the density difference that correspond to the pair 
amplitude nodes of the LO condensate. While the profile at t = 3 after a spin-flip excitation 
in figure d2l might seem similar, the peak of the density difference continues to move outward 
and the profile at t = 4 has three peaks (including the one at z - -0.45). While the local 
spin excitation develops into an oscillating density difference profile, the oscillation does not 
guarantee an LO condensate. 

Next we examine the effect of the location of the spin flip. In figure [Q} we plot the 
density profile after the spin flip at two different z c off the center, z c = -0.3125 and -0.625, 
respectively. The resulting density difference profile at t = does not have a peak at z c , 
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Figure 13. Density distributions for the ground state for (N-[,N]) = (9,7) atoms in the trap 
with L = 32, A/ J = 1 = 1024/L 2 and U/J = -8 = -256/L. 




Figure 14. Density distributions at different times after the spin-flip excitation at t — with 
z c = -0.3125 (left) and z c = -0.625 (right), plotted against the normalized coordinate for 
L = 32, w = 1/16, U/J = -8 and A/ J = 1. 



because in calculating the excited state (fTSt . the amplitude of the spin-flipped wave function 
S^l^o), < v F |5 / _ 5 / + | l r'o) J depends on the site number /. As the state evolves in time, the peak 
splits into two and they move to the left and right, as in the case of the spin flip at z = 0. 
However, the movement of the density difference peak to the direction of z > is considerably 
faster than that to the direction of z < 0, as the kinetic energy of the atoms is greater for z ~ 0. 

In figure El as the spin flip location is removed from the trap center, we observe that the 
change of the minority population in the side of the spin flip (z < 0) becomes less pronounced, 
and that in the other side of the trap, the densities of both components stay almost the same. 
This is because the spin-imbalanced gas in the z < side pushes the spin-balanced region 
(z > 0) away, thus compensating for the reduced kinetic energy by increasing the potential 
energy in the trap, but infiltrating this region only slowly. 

Finally we look at the effect of the interaction strength. For a smaller \U\, which 
corresponds to a smaller |g 1D | and a larger |aJ D |, an increase in energy due to the spin flip 
is smaller, but at the same time, the effect of Fermi pairing might become weaker. 

In figure [15] we show the evolution of the system for U/J = -2 and A/7 = 1, after the 
spin flip at the trap center (z c = 0). Observe that, while the atoms are extended over a wider 
region due to the weaker attraction, the speed of the movement of the two density-difference 
peaks is almost the same as in the case of U/J = -8. The insensitivity of the speed to the 
interaction strength, which is also observed for other spin flip parameters, suggests that the 
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Figure 15. Density distributions at different times after a spin-flip excitation at t — plotted 
against the normalized coordinate, for L = 32, z c = 0, w = 1/16, U/J — —2 and A// = 1. 

momenta of the excess up-spin atoms are determined primarily from the shape of the center 
peak whose width is independent of \U\, and not much affected by U, which determines 
the pairing strength. Indeed, if \U\ is kept constant and w is increased, the width of the 
initial distribution of the population imbalance increases, and the speed of the diffusion of the 
density difference rapidly decreases; this suggests that the speed may be determined from the 
uncertainty principle in the initial peak where the excess up-spin atoms are introduced by the 
spin flip. 

4. Conclusion 

We have studied the ground state and dynamics of harmonically confined ID Fermi gases 
with population imbalance. We have discretized the ID system and applied the density- 
matrix renormalization group method, which is numerically exact at least for the ground 
state in our setup, and also allows the simulation of the dynamics of large quantum systems 
at an unprecedented precision. For the ground state, we have demonstrated that for finite 
imbalance and attractive interaction, a Larkin-Ovchinnikov type condensate forms at the trap 
center, for a wide range of the polarization parameter up to almost unity. Qualitatively this 
result is consistent with both numerical studies O O 031 |20l EH E21 E3l, and studies that 
utilize the exact solution for a uniform system together with local density approximation 
fM 12211281 |29l El. 

Beyond a critical value of P, the LO condensate close to the trap center phase-separates 
from the spin-polarized Fermi gas at the periphery, regardless of the strength of interaction. 
By studying the dependence of the eigenvalues of two-body density matrix on P and N, we 
have shown that the oscillating pair correlation originates from quasi-condensation of atom 
pairs, to which almost all the minority atoms contribute. The FFLO-like feature is found to 
be rather insensitive to the asymmetry, which again supports the validity of the local density 
approximation. 

We have studied the dynamics after a spin-flip excitation by time-dependent DMRG 
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simulations. Of particular interest is whether the density difference oscillates as in the LO 
condensate in the ground state. We have found that, while the snapshots of the density 
difference show two or more peaks that originate from a single peak after the localized spin 
flip, the peaks move in time and are not related with an FFLO-like state. 
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